Reaction-diffusion system with self-organized critical behavior

We describe the construction of a conserved reaction-diffusion system that exhibits self-organized critical (avalanche-like) behavior under the action of a slow addition of particles. The model provides an illustration of the general mechanism to generate self-organized criticality in conserving systems. Extensive simulations in d = 2 and 3 reveal critical exponents compatible with the universality class of the stochastic Manna sandpile model. PACS. 05.70.Ln Nonequilibrium and irreversible thermodynamics – 05.65.+b Self-organized systems Since the introduction of the Bak, Tang, and Wiesenfeld sandpile model [1], the concept of self-organized criticality (SOC) [2] has witnessed a real explosion of activity, covering both the description of new models and the proposal of several theoretical approaches, aiming at an understanding of SOC phenomena in terms of standard statistical mechanical concepts. At this respect, it has been shown that SOC in sandpile models is related to the behavior of absorbing-state phase transitions (APT) with many absorbing states [3,4]. Indeed, this very idea is underlying a recipe proposed for the construction of SOC models [5]: any cellular automata, defined with conserved microscopic rules, and possessing many absorbing states, will display SOC behavior if slowly driven with the addition of energy/particles at a rate h and dissipation at a rate : i.e., in the double limit h → 0, → 0, with h/ → 0 [6]. This mechanism is easily seen at work in all standard sandpile models proposed so far [1,7,8]. Given that most SOC systems are defined in terms of sandpile-like models (with the exception of the forestfire [9] and extremal [10] models), it becomes all the most interesting to explore the possibility of applying the recipe of reference [5] to models of a qualitatively different sort. In this paper we consider a reaction-diffusion (RD) model showing an APT that conserves the total number of particles [11,12]. This model exhibits a non-equilibrium phase transition in the same universality class of fixed energy stochastic sandpiles [4,12]. Here, we show that implementing the slow driving condition, the model reaches a stationary state with an avalanche-like reaction activity with critical properties. By measuring usual magnitudes characterizing the SOC behavior, we compare the model with standard slowly driven sandpiles. The critical exponents measured confirm the shared universality class with a e-mail: [email protected] stochastic sandpiles, and provide a vivid illustration of the SOC generating mechanism [5]. We focus on the two species RD system [11,12], recently proposed to describe APT coupled to a nondiffusive conserved field [13]. The RD system is defined by the following set of reaction steps: B → A with rate k1, (1) B +A→ 2B with rate k2. (2) In this system, B particles diffuse with diffusion rate DB, and A particles do not diffuse, that is, DA = 0. From the rate equations (1) and (2), it is clear that the dynamics conserves the total number of particles N = NA + NB, where Ni is the number of particles i = A,B. In this model, the dynamics is exclusively due to B particles, that we identify as active particles. A particles do not diffuse and cannot generate spontaneously B particles. More specifically, A particles can only move via the motion of B particles that later on transform into A through equation (2). This implies that any configuration devoid of B particles is an absorbing state in which the system is trapped forever. The number of these absorbing states is infinite – in the thermodynamic limit – corresponding to all the possible redistributions of N particles of type A in the system. This RD process exhibits a phase transition from an active phase (with an everlasting activity of B particles) to an absorbing phase (no B particles) for a critical value ρ = ρc of the total particle density [12]. Here, we define a driven-dissipative version of the RD model by applying the recipe of references [4,5]. On hypercubic lattices of size L with open boundary conditions, each site i stores a number ai of A particles and bi of B particles. The occupation numbers ai and bi can have any integer value, including ai = bi = 0, that is, void sites with no particles. The model is thus representing the 584 The European Physical Journal B dynamics of bosonic particles. The initial configuration is constructed by randomly distributing a number N0 of A particles in the lattice. The initial occupation numbers ai have a Poissonian distribution, while bi = 0,∀i. Any configuration is stable whenever it fulfills this condition, i.e., in the absence of B particles. The system is driven by adding one B particle to a randomly chosen site. A state with at least one B particle is called active. Active states evolve in time according to the following update rules, that mimic the diffusion and reaction steps in the RD system: I) Diffusion: on each lattice site, each B particle moves into a randomly chosen nearest neighbor site with probability 2d/(2d + 1), and stays in the same site with probability 1/(2d+1); this results in an effective diffusivity DB = 1/(2d+ 1). II) After all sites have been updated for diffusion, we perform the reactions: a) On each lattice site, each B particle is turned into an A particle with probability r1. b) At the same time, each A particle becomes a B particle with probability 1−(1−r2)i , where bi is the total number of B particles in that site. This corresponds to the average probability for an A particle of being involved in the reaction (2) with any of the B particles present on the same site. The probabilities r1 and r2 are related to the reaction rates k1 and k2 defined in equations (1) and (2). In general, we have that ri(ki = 0) = 0, ri(ki = ∞) = 1, and ri is an increasing function of ki. The analytic expression of ri as a function of ki is presumably quite complex and nontrivial. However, as we will argue later, the knowledge of the precise relationship between ri and ki is not necessary, since the critical behavior of the model should be independent of the exact values of the parameters ri selected. B particles on boundary sites may choose to diffuse out of the lattice. In this case, the particle is removed out of the system, contributing to the dissipation. The system is updated in parallel until there are no more B particles and it is again in an absorbing state. During the dynamic evolution, the addition of new B particles is suspended; this of course corresponds to the slow-driving condition. The sequence of updates in the system (from the time we introduce a new B particle until a stable state is reached) is interpreted as an avalanche. We characterize avalanches by their size s and their duration t. The size of an avalanche is defined as the number of B particles present in the system at each time step, summed over all the parallel updates required to reach a new stable state. The duration of an avalanche is defined as the total number of parallel updates performed during the avalanche. In the slow driving perspective, the existence of a critical stationary state is easily understood. Particles are added only in the absence of activity (ρ < ρc), while dissipation acts only during activity (ρ > ρc). This implies that ∂tρ always drives the system toward ρc, that in the thermodynamic limit is the only possible stationary value of the density [4,5]. We have performed numerical simulations of this model in dimensions 2 and 3, with system sizes ranging from L = 64 to L = 1024 in d = 2, and from L = 74 to L = 280 in d = 3. The reaction rates ri reported here are r1 = 0.3 and r2 = 0.4 in d = 2, and r1 = 0.4 and r2 = 0.5 0 20 40 60 80 100 10 2 10 3 10 4 10 5 10 0 10 1 10 2 10 3 10 4 10 5